R Notebook

This document presents the subset of the figures used for paper about monitoring.

Load packages

library(MASS)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggmap)
library(scatterpie)
library(rgdal)
library(multcompView)
library(car)
library(ggpmisc)
library(chisq.posthoc.test)
library(vcd)
library(grid)
library(ggpattern)
library(car)
library(gridExtra)
library(cowplot)

Sampling data correspond to the data collected with kobo and previously cleaned with the script 1_preprocesamiento_datos_kobo.Rmd.

Load data

# load data
muestreo_tidy<-read.delim("../data/kobo/muestreo_dic2020_tidy.txt", header = TRUE)

#Add column block to muestreo_tidy
muestreo_tidy<-muestreo_tidy%>% mutate(block = case_when(plot %in% c(1,2,3) ~1,
                 plot %in% c(1,2,3) ~1,
                 plot %in% c(4,5,6) ~2,
                 plot %in% c(7,8,9) ~3,
                 plot %in% c(10,11,12) ~4,
                 plot %in% c(13,14,15) ~5,
                 plot %in% c(16,17,18) ~6,
                 plot %in% c(19,20,21) ~7,
                 plot %in% c(22,23,24) ~8,
                 plot %in% c(25,26,27) ~9,
                 plot %in% c(28,29,30) ~10,
                 plot %in% c(31,32,33) ~11,
                 plot %in% c(34,35,36) ~12,
                 plot %in% c(37,38,39) ~13,
                 plot %in% c(40,41,42) ~14,
                 plot %in% c(43,44,45) ~15,
                 plot %in% c(46,47,48) ~16)) %>%
                mutate(block=as.factor(block)) # store as factor, instead of numeric (numbers are IDs)

parcelas_tidy<-read.delim("../data/kobo/parcelas_dic2020_tidy.txt", header = TRUE)

# pivot long plots data to have health data as a single variable
parcelas_long<-pivot_longer(parcelas_tidy, 
                            cols = healthy:worm, 
                            names_to = "tree_health_simplified",
                            values_to = "n_trees")

Data analyzed here correspond only to the trees that were approved during the validation by manually reviewing the photographs in kobotoolbox. Total of 1778 trees sampled, 1765 were approved in the validation.

muestreo_tidy<- filter(muestreo_tidy, X_validation_status=="validation_status_approved")

# write.csv(muestreo_tidy,  file="../../../../Desktop/muestreo_tidy.csv")

Color palettes:

# Make a nice color pallete and legend order for all plots

my_cols=c("darkgreen", 
              "darkred", 
              "orangered1", 
              "cadetblue", 
              "tan", 
              "beige", 
            #  "burlywood4", 
              "coral", 
              "aquamarine3", 
              "gray70", 
              "black")

desired_order=c("healthy", 
                "ozone", 
                "ozone_and_other", 
                "others_combined", 
                "drougth", 
                "fungi", 
             #   "insect", 
                "worm", 
                "acid_rain", 
                "other", 
                "dead")

desired_names=c("healthy", 
                "ozone", 
                "ozone and other", 
                "others combined", 
                "drougth", 
                "fungi", 
             #   "insect", 
                "worm", 
                "acid rain", 
                "other", 
                "dead")

spanish_labels=c("Sano", 
                  "Ozono", 
                  "Ozono y otros", 
                  "Otros combinados no-ozono", 
                  "Sequía", 
                  "Hongos", 
               #   "Insectos", 
                  "Gusano de seda", 
                  "Lluvia acida", 
                  "Otro", 
                  "Muerto")

# For ozone damage percentage 
 my_cols2<-c("darkgreen", "gold2", "chocolate1", "orangered", "red4", "darkorchid4")
 
desired_order_percentage<-c("0%","less than 10%", "10 to 40%", "40 to 50%", "50 to 70%", "more than 70%")

Multiplot fun:

# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

Configure google api for maps:

# code adapted from https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-g0ogle-map.html

## configure google api

# You first need to register your api key in https://cloud.google.com/maps-platform/#get-started and follow instructions. The geocoding API is a free service, but you nevertheless need to associate a credit card with the account. Please note that the Google Maps API is not a free service. There is a free allowance of 40,000 calls to the geocoding API per month, and beyond that calls are $0.005 each.
# after you obtain your api, save it in /scripts/api_key.api (not shown in this repo por obvious reasons).

# if you get the following error when running get_map():

#"Error in aperm.default(map, c(2, 1, 3)) : 
#  invalid first argument, must be an array " 

# check this troubleshooting: https://rgraphgallery.blogspot.com/2013/04/rg-plot-pie-over-google-map.html

##  load and register api
api <- readLines("api_key.api")
register_google(key = api)

Map and monitoring figures presented in the paper:

Figure 2

Plot 2a: PNDL location on CDMX map

# get cdmx shape
CDMX<-readOGR(dsn="../data/spatial", layer="CDMX")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/veronicareyesgalindo/Documents/GitHub/monitoreo-oyameles/data/spatial", layer: "CDMX"
## with 1 features
## It has 8 fields

CDMX<-fortify(CDMX)

# get PNDL shape
PNDL<-readOGR(dsn="../data/spatial", layer="Desierto_Leones_Geo_ITRF08")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/veronicareyesgalindo/Documents/GitHub/monitoreo-oyameles/data/spatial", layer: "Desierto_Leones_Geo_ITRF08"
## with 1 features
## It has 14 fields

PNDL<-fortify(PNDL)

# get background map
sat_map = get_map(location = c(lon = -99.133549, lat = 19.3), zoom = 10, maptype = 'terrain-background', source = "google")

## plot
p_a<-ggmap(sat_map) + 
            geom_polygon(data = CDMX,
                         aes(x = long, y = lat, group = group),
                         color="black", fill=NA, size=1.5) +
            geom_polygon(data = PNDL,
                         aes(x = long, y = lat, group = group),
                         color="red", fill=NA, size=1.5) +
            geom_point(aes(x=-98.95, y=19.6), 
                       shape=0, stroke=2, size=5, color="black") +
            geom_point(aes(x=-98.95, y=19.55), 
                       shape=0, stroke=2, size=5, color="red") +
            geom_text(aes(label="CDMX", x=-98.87, y=19.6), 
                      color="Black", fontface="bold", size=5) +
            geom_text(aes(label="PNDL", x=-98.87, y=19.55), 
                      color="Black", fontface="bold", size=5) +
            theme(text = element_text(size = 20))+
  ggtitle("a)")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Plot 2b: Satellite image and surroundings of the PNDL

# get background map
sat_map = get_map(location = c(lon = -99.30, lat = 19.31), zoom = 13, maptype = 'satellite', source = "google")

## add towns names
towns<-data.frame(nombre=c("San Bartolo Ameyalco", 
                           "Santa Rosa Xochiac", 
                           "San Mateo Tlaltenango"),
                  long=c(-99.270, -99.29, -99.276),
                  lat=c(19.333, 19.325, 19.346))



## plot
p_b<-ggmap(sat_map) + 
            geom_polygon(data = PNDL,
                         aes(x = long, y = lat, group = group),
                         color="red", fill=NA, size=1.5) +
            geom_point(data=towns, aes(x=long, y=lat), colour="red", size=1.5) +
            geom_text(data=towns, aes(label=nombre, x=long, y=lat), 
                      color="white", fontface="bold",
                      size=5, nudge_y=0.003) +
  # add Cruz de Coloxtitla (CX), and Convento (Cn) landmarks
            geom_text(aes(label="X", x=-99.3014, y=19.286068), 
                      color="white", fontface="bold", size=4) +
            geom_text(aes(label="C", x=-99.31, y=19.3133), 
                      color="white", fontface="bold", size=4) +
            theme(text = element_text(size = 20))+
  ggtitle("b)")

Plot 2c: This is the distribution of the 48 plots:

## plot map
# get map
sat_map = get_map(location = c(lon = -99.3060, lat = 19.2909), zoom = 14, maptype = 'satellite', source = "google")

# plot sampled plots
p_c <-  ggmap(sat_map)
p_c <- p_c + geom_point(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude),
                      color="red") +
          geom_text(data=parcelas_tidy,
                      aes(x=X_coordinates_longitude,
                          y=X_coordinates_latitude,
                          label=plot),
                      color="white",
                     check_overlap = TRUE,
                      hjust = 0, vjust=1, nudge_x = 0.0005,
                 size= 5) +
    theme(text = element_text(size = 20))+
  ggtitle("c)")
p_c 

Plot 2d: Distribution of tree health status by plot

The following figure shows the total number of trees sampled in each 10x10 m plot, and how many of these are under some category of damage:

parcelas_HOO<-parcelas_long%>%
  filter(tree_health_simplified ==  "ozone" )
  sum(parcelas_HOO$n_trees)

## [1] 125

parcelas_HOO<-parcelas_long%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" | tree_health_simplified == "ozone_and_other" )


p_d <- ggplot(parcelas_HOO, aes(x=plot, y=n_trees,     fill=tree_health_simplified)) +
  geom_bar(stat="identity") +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud")
  
p_d <- p_d + theme_bw() +
  labs(x="Parcelas", y= "Número de árboles") +
  theme(text = element_text(size = 20)) +
  ggtitle("d)")+
  coord_flip()

Multiplot Fig 2

fig_2<-multiplot(p_a, p_c, p_b, p_d, cols=2)

ggsave("../outputs/Fig_2.png",
       dpi = 600)

ANOVAs related to Figure 3

One-way ANOVAS and 2-way ANOVA

#### subset of data only with reforestation, health and exposition
refo_expo<-muestreo_tidy %>% 
              dplyr::select(plot, block, tree_health_simplified, reforested, tree_exposition) %>% 
           # discard other typos of damage
             filter(tree_health_simplified == "healthy" |
                    tree_health_simplified == "ozone" |
                    tree_health_simplified == "ozone_and_other") %>% 
          # make "ozone" and "ozone and other" a single category = damaged
             mutate(tree_health_simplified = replace(tree_health_simplified, tree_health_simplified == "ozone", "damaged")) %>% 
             mutate(tree_health_simplified = replace(tree_health_simplified, tree_health_simplified == "ozone_and_other", "damaged"))
      


### ONE-way ANOVA
library(car)

#### By reforestation
summary_by_ref<- refo_expo %>%  
  group_by(plot, block, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

## `summarise()` has grouped output by 'plot', 'block', 'reforested'. You can
## override using the `.groups` argument.

#now we have n for each category
head(summary_by_ref)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_ref<-transform(summary_by_ref, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_ref) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_ref<-summary_by_ref %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore

# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK





## By exposition
summary_by_exp<- refo_expo %>%  
  group_by(plot, block, tree_exposition,  tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

## `summarise()` has grouped output by 'plot', 'block', 'tree_exposition'. You can
## override using the `.groups` argument.

#now we have n for each category
head(summary_by_exp)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_exp<-transform(summary_by_exp, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_exp) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_exp<-summary_by_exp %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore



#### One way ANOVAS

# for reforestation
Anova(lm(summary_by_ref$perc.damage.log ~ summary_by_ref$reforested + summary_by_ref$block))

# for exposition
Anova(lm(summary_by_exp$perc.damage.log ~ summary_by_exp$tree_exposition + summary_by_exp$block))

#### For TWO-way ANOVA

### Summarise and estimate n and % of DAMAGED trees (normalized) by category and plot

summary_by_both<- refo_expo %>%  
  group_by(plot, block, tree_exposition, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

## `summarise()` has grouped output by 'plot', 'block', 'tree_exposition',
## 'reforested'. You can override using the `.groups` argument.

#now we have n for each category
head(summary_by_both)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_both<-transform(summary_by_both, 
                           perc.damage = ave(n,
                                      plot,
                                      FUN = prop.table))  %>% 
                      
 # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_both) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_both<-summary_by_both %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore


# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK

head(summary_by_both) 

##### TWO WAY ANOVA By BOTH exposition and tree origin 

### get data in correct format
summary_by_both<- refo_expo %>%  
  group_by(plot, block, tree_exposition, reforested, tree_health_simplified) %>% 
  # count how many individuals there are in each category
  summarise(n = n())

## `summarise()` has grouped output by 'plot', 'block', 'tree_exposition',
## 'reforested'. You can override using the `.groups` argument.

#now we have n for each category
head(summary_by_both)

# Calculate percentage by group and plot ((summing all rows of each plot = 1)
summary_by_both<-transform(summary_by_both, 
                           perc.damage = ave(n,
                                             plot,
                                             FUN = prop.table))  %>% 
  
  # normalize it
  mutate(perc.damage.log = log(perc.damage))

# the percentage of healthy and damaged trees for each exposed covered/exposed category is shown. 
head(summary_by_both) 

# Next we need to remove the "healthy rows" (we only need the % of damage, not the % of healthy)
summary_by_both<-summary_by_both %>%
  filter(tree_health_simplified == "damaged") %>%
  select(-tree_health_simplified) # this column is not needed anymore


# The percentages now don't sum 1 because we are only looking at the % of damaged trees by category, NOT to the % of healthy trees. This is OK

head(summary_by_both) 

Figure S10

### TWO-way ANOVA
Anova(lm(summary_by_both$perc.damage.log ~ summary_by_both$tree_exposition * summary_by_both$reforested + summary_by_both$block))

Figure 3 Boxplot of data used for ANOVA

# Boxplot tree expostion and reforestation origin on percentage of damage
p3<- ggplot(summary_by_both, aes(x = factor(tree_exposition), 
                                 y = perc.damage.log, color = factor(reforested))) +
  geom_boxplot() +
  ylab("Log-transformado del porcentaje de daño por árbol") +
  scale_color_manual(name= "Origen del árbol",
                     breaks = c("no", "yes"),
                     labels = c("Regeneración natural", "Reforestedo"),
                     values=c("#9f964b", "#9a5ea1")) + # nice colour blind colors
  
  scale_x_discrete(name="Exposición del árbol",
                   breaks = c("cover", "exposed"),labels = c("Con nodriza", "Expuesto")) + 
  theme_bw()+ 
  theme(text = element_text(size = 20))

p3

Figure 4

Plot 4a

condition_HOO<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |
           tree_health_simplified == "ozone_and_other" )


p <- filter(condition_HOO, tree_heigth<15, tree_nodes>0) %>% 
     ggplot(.) +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") 
theme_bw()

## List of 97
##  $ line                      :List of 6
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ lineend      : chr "butt"
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ rect                      :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : chr "black"
##   ..$ linewidth    : num 0.5
##   ..$ linetype     : num 1
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ text                      :List of 11
##   ..$ family       : chr ""
##   ..$ face         : chr "plain"
##   ..$ colour       : chr "black"
##   ..$ size         : num 11
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 0
##   ..$ lineheight   : num 0.9
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ title                     : NULL
##  $ aspect.ratio              : NULL
##  $ axis.title                : NULL
##  $ axis.title.x              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.75points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.top          :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.75points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.bottom       : NULL
##  $ axis.title.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.75points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y.left         : NULL
##  $ axis.title.y.right        :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.75points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text                 :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey30"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.2points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.top           :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 2.2points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.bottom        : NULL
##  $ axis.text.y               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 1
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 2.2points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y.left          : NULL
##  $ axis.text.y.right         :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 0points 2.2points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks                :List of 6
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ axis.ticks.x              : NULL
##  $ axis.ticks.x.top          : NULL
##  $ axis.ticks.x.bottom       : NULL
##  $ axis.ticks.y              : NULL
##  $ axis.ticks.y.left         : NULL
##  $ axis.ticks.y.right        : NULL
##  $ axis.ticks.length         : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ axis.ticks.length.x       : NULL
##  $ axis.ticks.length.x.top   : NULL
##  $ axis.ticks.length.x.bottom: NULL
##  $ axis.ticks.length.y       : NULL
##  $ axis.ticks.length.y.left  : NULL
##  $ axis.ticks.length.y.right : NULL
##  $ axis.line                 : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x               : NULL
##  $ axis.line.x.top           : NULL
##  $ axis.line.x.bottom        : NULL
##  $ axis.line.y               : NULL
##  $ axis.line.y.left          : NULL
##  $ axis.line.y.right         : NULL
##  $ legend.background         :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.margin             : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing            : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ legend.spacing.x          : NULL
##  $ legend.spacing.y          : NULL
##  $ legend.key                :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.key.size           : 'simpleUnit' num 1.2lines
##   ..- attr(*, "unit")= int 3
##  $ legend.key.height         : NULL
##  $ legend.key.width          : NULL
##  $ legend.text               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.text.align         : NULL
##  $ legend.title              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.title.align        : NULL
##  $ legend.position           : chr "right"
##  $ legend.direction          : NULL
##  $ legend.justification      : chr "center"
##  $ legend.box                : NULL
##  $ legend.box.just           : NULL
##  $ legend.box.margin         : 'margin' num [1:4] 0cm 0cm 0cm 0cm
##   ..- attr(*, "unit")= int 1
##  $ legend.box.background     : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ legend.box.spacing        : 'simpleUnit' num 11points
##   ..- attr(*, "unit")= int 8
##  $ panel.background          :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.border              :List of 5
##   ..$ fill         : logi NA
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.spacing             : 'simpleUnit' num 5.5points
##   ..- attr(*, "unit")= int 8
##  $ panel.spacing.x           : NULL
##  $ panel.spacing.y           : NULL
##  $ panel.grid                :List of 6
##   ..$ colour       : chr "grey92"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.grid.major          : NULL
##  $ panel.grid.minor          :List of 6
##   ..$ colour       : NULL
##   ..$ linewidth    : 'rel' num 0.5
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.grid.major.x        : NULL
##  $ panel.grid.major.y        : NULL
##  $ panel.grid.minor.x        : NULL
##  $ panel.grid.minor.y        : NULL
##  $ panel.ontop               : logi FALSE
##  $ plot.background           :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : chr "white"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ plot.title                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 5.5points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.title.position       : chr "panel"
##  $ plot.subtitle             :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0points 0points 5.5points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.caption              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : num 1
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 5.5points 0points 0points 0points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.caption.position     : chr "panel"
##  $ plot.tag                  :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.tag.position         : chr "topleft"
##  $ plot.margin               : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
##   ..- attr(*, "unit")= int 8
##  $ strip.background          :List of 5
##   ..$ fill         : chr "grey85"
##   ..$ colour       : chr "grey20"
##   ..$ linewidth    : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ strip.background.x        : NULL
##  $ strip.background.y        : NULL
##  $ strip.clip                : chr "inherit"
##  $ strip.placement           : chr "inside"
##  $ strip.text                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey10"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
##   .. ..- attr(*, "unit")= int 8
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.x              : NULL
##  $ strip.text.x.bottom       : NULL
##  $ strip.text.x.top          : NULL
##  $ strip.text.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.y.left         :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.y.right        : NULL
##  $ strip.switch.pad.grid     : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  $ strip.switch.pad.wrap     : 'simpleUnit' num 2.75points
##   ..- attr(*, "unit")= int 8
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi TRUE
##  - attr(*, "validate")= logi TRUE

p4_a <- p + geom_histogram(aes(x=tree_nodes, fill=tree_health_simplified))  +
    labs(x="Edad del árbol (años)", y= "Número de árboles") +
    theme(text = element_text(size = 20)) +
    theme(plot.title = element_text(lineheight=1.1))+
  ggtitle("a)")
p4_a

Plot 4c

## base data
# Definir plantas sanas y dañadas por otra cosa que no fuera ozono
# cond_PO<- se  refiere a condition Percentage damage by Ozone 
cond_PO<-as_data_frame(muestreo_tidy)

## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` (with slightly different semantics) to convert to a
##   tibble, or `as.data.frame()` to convert to a data frame.

# Asignar 0% de daño por ozono a los árboles healthy
cond_PO$ozone_damage_percentage = ifelse(cond_PO$tree_health == "healthy", "0%", cond_PO$ozone_damage_percentage)

# Filtrar por porcentaje de daño
condition_PO<-cond_PO%>%
  filter(ozone_damage_percentage == "0%" | ozone_damage_percentage == "less than 10%" | ozone_damage_percentage == "10 to 40%" | ozone_damage_percentage == "40 to 50%"| ozone_damage_percentage == "50 to 70%" | ozone_damage_percentage == "more than 70%")

condition_PO$ozone_damage_percentage <- as.factor(condition_PO$ozone_damage_percentage)

library(sf)
library(ggpattern)
library(magrittr)
library(ggplot2)
# Plot
p_od<- condition_PO %>% filter(!is.na(ozone_damage_percentage)) %>%
  ggplot(., aes(x = tree_nodes, fill = ozone_damage_percentage)) +
  geom_bar_pattern(aes(pattern = ozone_damage_percentage),
                   color = "black",
                   pattern_fill = "black",
                   pattern_density = 0.1,
                   pattern_spacing = 0.025,
                   pattern_key_scale_factor = 0.6) +
  scale_fill_manual(values= my_cols2, 
                    breaks = desired_order_percentage,
                    labels = c("0%","menos de 10%", "10 a 40%", "40 a 50%",
                               "50 a 70%", "más de 70%"),
                    name= "Porcentaje del árbol\n dañado por ozono") +
  scale_pattern_manual(name= "Porcentaje del árbol\n dañado por ozono",
                       values =c("none","none", "stripe", "none",
                                 "none", "none"),
                       labels = c("0%","menos de 10%", "10 a 40%", "40 a 50%",
                               "50 a 70%", "más de 70%"))+
  theme_bw()+ theme(text = element_text(size = 20)) 


p_od

## Warning: Removed 147 rows containing non-finite values (`stat_count()`).

p4_c <- p_od +
  labs(x="Edad del árbol (años)", y= "Número de árboles")  +
  theme(legend.title.align = 0.5)+
  theme(text = element_text(size = 20)) +
  theme(plot.title = element_text(lineheight=1.1))+
  ggtitle("c)")

p4_c

## Warning: Removed 147 rows containing non-finite values (`stat_count()`).

Plot 4b

# Filtrar por categoría de daño
condition_HOO<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" | tree_health_simplified == "ozone_and_other" )

condition_HOO$tree_health_simplified <- as.factor(condition_HOO$tree_health_simplified)



# Data distribution

# Los datos tienen a graficar es el número de nodos para cada categoria de salud.
# Los datos son continuos discretos, por lo tanto el analisis a seguir para buscar diferencias entre los grupos son:

shapiro.test(condition_HOO$tree_nodes) #NO HAY NORMALIDAD 

## 
##  Shapiro-Wilk normality test
## 
## data:  condition_HOO$tree_nodes
## W = 0.94294, p-value < 2.2e-16

# Procedo a hacer una prueba no paramétrica (kruskal)

# Debe tener Homocedasticidad: la varianza debe de ser constante entre todos los grupos.
# If the p-value is less than our chosen significance level, we can reject the null hypothesis and conclude that we have enough evidence to state that the variance among the groups is not equal.
leveneTest(sqrt(tree_nodes) ~ tree_health_simplified, data = condition_HOO) #Opera con medias

# As the p value obtained from the Shapiro-Wilk test and Levene’s test is significant (p < 0.05), we conclude that the data is not normally distributed and does not have equal variance. Kruskal-Wallis test is more appropriate for analyzing differences.

kruskal.test(sqrt(tree_nodes) ~ tree_health_simplified, data = condition_HOO) #opera con mediana

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sqrt(tree_nodes) by tree_health_simplified
## Kruskal-Wallis chi-squared = 265.57, df = 2, p-value < 2.2e-16

#As the p value obtained from the Kruskal-Wallis test test is significant , we conclude that there are significant differences.

# For the Kruskal-Wallis test, epsilon-squared is a method of choice for effect size measurement. The epsilon-squared is 0.74 and suggests a very strong effect of plant varieties on yield

# calculate effect size. Hay efecto relativamente fuerte
library(rcompanion)
epsilonSquared(x = sqrt(condition_HOO$tree_nodes), g = condition_HOO$tree_health_simplified)

## epsilon.squared 
##           0.239

# https://peterstatistics.com/CrashCourse/3-TwoVarUnpair/NomOrd/NomOrd3c.html
# To know which plant varieties are significantly different from each other, we will perform the Dunn’s test as post-hoc test for significant Kruskal-Wallis test. As there are multiple comparisons, we will correct the p values using Benjamini-Hochberg FDR method for multiple hypothesis testing at a 5% cut-off. The other tests that can be used for post-hoc test includes Conover and Nemenyi tests. Dunn’s test is more appropriate post-hoc than the Mann-Whitney U test for significant Kruskal-Wallis test as it retains the rank sums of the Kruskal-Wallis.

#poshoc para saber que grupos difieren
pairwise.wilcox.test(x = sqrt(condition_HOO$tree_nodes), g = condition_HOO$tree_health_simplified, p.adjust.method = "holm" )

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  sqrt(condition_HOO$tree_nodes) and condition_HOO$tree_health_simplified 
## 
##                 healthy ozone
## ozone           <2e-16  -    
## ozone_and_other <2e-16  0.012
## 
## P value adjustment method: holm

# Perorm pairwise comparisons
library(ggpubr)
my_comparisons <- list( c("healthy", "ozone"), c("healthy", "ozone_and_other"), c("ozone", "ozone_and_other") )


p4_b<-condition_HOO%>%
   ggplot(aes(y= tree_health_simplified, x= tree_nodes))+
        scale_color_manual(values=  my_cols, 
                           labels= spanish_labels,
                    name= "")+
        geom_point(position="jitter",aes(color = tree_health_simplified), alpha=0.5, size= 2.5)+
  geom_boxplot(color="black", notch = F, alpha = 0.1)+
        xlab("Edad del árbol (años)")+ ylab("Número de árboles")+
  theme_bw()+
  ggtitle("b)")+
  theme(text = element_text(size = 20), axis.text.y=element_blank())+
  theme(plot.title = element_text(lineheight=1.1))+
  scale_x_continuous(breaks=c(5,10,15, 20, 25))+
  stat_compare_means(comparisons = my_comparisons, label = "p.signif")

p4_b

## Warning: Removed 147 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 147 rows containing non-finite values (`stat_signif()`).

## Warning: Removed 147 rows containing missing values (`geom_point()`).

Plot 4d

# Perorm pairwise comparisons
# Run Shapiro-Wilk test
shapiro.test(condition_PO$tree_nodes) #NO HAY NORMALIDAD

## 
##  Shapiro-Wilk normality test
## 
## data:  condition_PO$tree_nodes
## W = 0.94294, p-value < 2.2e-16

# Procedo a hacer un kruskal
# Debe tener homogeneidad. Si es mayor a 0.05 No hay evidencias en contra de la homogeneidad de varianzas. 
leveneTest(sqrt(tree_nodes) ~ ozone_damage_percentage, data = condition_PO, center = "median")

kruskal.test(sqrt(tree_nodes) ~ ozone_damage_percentage, data = condition_PO)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sqrt(tree_nodes) by ozone_damage_percentage
## Kruskal-Wallis chi-squared = 283.95, df = 5, p-value < 2.2e-16

epsilonSquared(x = sqrt(condition_PO$tree_nodes), g = condition_PO$ozone_damage_percentage)

## epsilon.squared 
##           0.255

#poshoc que grupos difieren

#order levels 
condition_PO$ozone_damage_percentage<- ordered(condition_PO$ozone_damage_percentage, levels=c("0%","less than 10%", "10 to 40%", "40 to 50%", "50 to 70%","more than 70%"))

#Prueba
pairwise.wilcox.test(x = sqrt(condition_PO$tree_nodes), g = condition_PO$ozone_damage_percentage, p.adjust.method = "bonferroni" )

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  sqrt(condition_PO$tree_nodes) and condition_PO$ozone_damage_percentage 
## 
##               0%      less than 10% 10 to 40% 40 to 50% 50 to 70%
## less than 10% < 2e-16 -             -         -         -        
## 10 to 40%     < 2e-16 0.6498        -         -         -        
## 40 to 50%     < 2e-16 0.2086        1.0000    -         -        
## 50 to 70%     < 2e-16 2.9e-06       0.0026    0.5273    -        
## more than 70% 1.7e-05 1.0000        1.0000    1.0000    0.2060   
## 
## P value adjustment method: bonferroni

my_comparisons <- list(c("0%", "less than 10%"), c("0%", "10 to 40%"), c("0%", "40 to 50%"),
                        c("0%", "50 to 70%"), c("0%", "more than 70%"),
                        c("less than 10%", "50 to 70%"), 
                       c("10 to 40%", "50 to 70%"))

# Plot
p4_d<-condition_PO%>% filter(!is.na(ozone_damage_percentage)) %>% 
   ggplot(aes(y=ozone_damage_percentage , x= tree_nodes))+
        scale_color_manual(values=  my_cols2,labels = c("0%","menos de 10%", "10 a 40%", "40 a 50%","50 a 70%", "más de 70%"))+
        geom_point(position="jitter",aes(color = ozone_damage_percentage), alpha=0.5, size= 2.5)+
        xlab("Edad del árbol (años)")+ ylab("Número de árboles")+
  labs(color = "")+
          geom_boxplot(color="black", notch = F, alpha = 0.1)+
  theme_bw()+
  ggtitle("d)")+
  theme(legend.title.align = 0.5)+
  theme(text = element_text(size = 20),axis.text.y=element_blank())+
  theme(plot.title = element_text(lineheight=1.1))+
  scale_x_continuous(breaks=c(5,10,15, 20, 25))+
  ggpubr::stat_compare_means(comparisons = my_comparisons, label="p.signif")


p4_d

## Warning: Removed 147 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 147 rows containing non-finite values (`stat_signif()`).

## Warning: Removed 147 rows containing missing values (`geom_point()`).

Multiplot Fig 4

multiplot(p4_a, p4_c, p4_b, p4_d, cols=2)

## Warning: Removed 147 rows containing non-finite values (`stat_count()`).

## Warning: Removed 147 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 147 rows containing non-finite values (`stat_signif()`).

## Warning: Removed 147 rows containing missing values (`geom_point()`).

## Warning: Removed 147 rows containing non-finite values (`stat_boxplot()`).

## Warning: Removed 147 rows containing non-finite values (`stat_signif()`).

## Warning: Removed 147 rows containing missing values (`geom_point()`).

Figure 5

Plot model all data WITHOUT cohort

condition_model_all<-dplyr::select(muestreo_tidy, contains(c("tree_heigth","tree_nodes","tree_health_simplified", "reforested", "tree_exposition", "plot")))%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |tree_health_simplified == "ozone_and_other" )%>%
  filter(tree_nodes < 15)

# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5all<-g0+geom_point(aes(colour=tree_health_simplified, 
                           shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified), fullrange =T)+
  scale_color_manual(values=  my_cols,
                     labels= spanish_labels,
                    name= "Estado de salud")+ 
  scale_shape_manual(name= "Estado de salud", labels= spanish_labels, values = c(15,16,17)) + 
  labs(x="Edad del árbol (años)", y= "Altura de los árboles")+
  theme_bw()+ 
  ggtitle("a)")
fig5all

Select cohort data

cohort_data<- dplyr::select(condition_HOO, contains(c("tree_health_simplified", "tree_nodes")))%>% filter(tree_health_simplified == "ozone")
summary(cohort_data)

##      tree_health_simplified   tree_nodes    
##  healthy        :  0        Min.   : 3.000  
##  ozone          :125        1st Qu.: 7.000  
##  ozone_and_other:  0        Median : 8.000  
##                             Mean   : 9.038  
##                             3rd Qu.:11.000  
##                             Max.   :19.000  
##                             NA's   :20

# We choose 7 & 11

Model all data with cohort 7-11

condition_model_all_C<-dplyr::select(muestreo_tidy, contains(c("tree_heigth","tree_nodes","tree_health_simplified", "reforested", "tree_exposition", "plot")))%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone" |tree_health_simplified == "ozone_and_other" )%>%
  filter(tree_nodes %in% (7:11))%>%
  filter(tree_nodes < 15)


########################
library(lme4)

########################

# Model all independent variables WITH health status interaction 
glm_1_coh<-glm(tree_heigth ~ tree_nodes, data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_1_coh)

## 
## Call:
## glm(formula = tree_heigth ~ tree_nodes, family = Gamma(link = "log"), 
##     data = condition_model_all_C)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6684  -0.6477  -0.1315   0.2175   1.9563  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.82092    0.24866  -3.301  0.00105 ** 
## tree_nodes   0.22820    0.02815   8.108 6.36e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.5628858)
## 
##     Null deviance: 235.07  on 401  degrees of freedom
## Residual deviance: 200.90  on 400  degrees of freedom
## AIC: 1642.7
## 
## Number of Fisher Scoring iterations: 5

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_1_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested data
glm_2_coh<-glmer(tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_2_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition +  
##     reforested + (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1525.5   1557.5   -754.7   1509.5      394 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5686 -0.7133 -0.0860  0.4988  4.9872 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.09391  0.3064  
##  Residual             0.33133  0.5756  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)                           -1.02639    0.24521  -4.186 2.84e-05 ***
## tree_nodes                             0.22537    0.02568   8.776  < 2e-16 ***
## tree_health_simplifiedozone            0.22387    0.09537   2.348   0.0189 *  
## tree_health_simplifiedozone_and_other  0.21017    0.07834   2.683   0.0073 ** 
## tree_expositionexposed                 0.11226    0.06787   1.654   0.0981 .  
## reforestedyes                         -0.32335    0.10358  -3.122   0.0018 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____ tr_xps
## tree_nodes  -0.913                            
## tr_hlth_smp -0.054 -0.071                     
## tr_hlth_s__ -0.195 -0.006  0.468              
## tr_xpstnxps -0.035 -0.063 -0.093  0.071       
## reforestdys -0.184  0.064 -0.073  0.021  0.015

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_2_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT tree exposition data
glm_3_coh<-glmer(tree_heigth ~ tree_nodes + tree_health_simplified + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00391099 (tol = 0.002, component 1)

summary(glm_3_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes + tree_health_simplified + reforested +  
##     (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1526.3   1554.2   -756.1   1512.3      395 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5830 -0.6956 -0.0782  0.5100  4.8349 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.09336  0.3056  
##  Residual             0.32926  0.5738  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)                           -1.01123    0.24581  -4.114 3.89e-05 ***
## tree_nodes                             0.22811    0.02572   8.867  < 2e-16 ***
## tree_health_simplifiedozone            0.23865    0.09525   2.505  0.01223 *  
## tree_health_simplifiedozone_and_other  0.20066    0.07833   2.562  0.01041 *  
## reforestedyes                         -0.32621    0.10383  -3.142  0.00168 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____
## tree_nodes  -0.918                     
## tr_hlth_smp -0.055 -0.080              
## tr_hlth_s__ -0.191 -0.002  0.479       
## reforestdys -0.183  0.065 -0.072  0.020
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00391099 (tol = 0.002, component 1)

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_3_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested and tree exposition data
glm_4_coh<-glmer(tree_heigth ~ tree_nodes+tree_health_simplified + tree_exposition + (1|plot) , data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_4_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition +  
##     (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1532.8   1560.8   -759.4   1518.8      395 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5286 -0.6923 -0.0918  0.5039  4.9279 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.1090   0.3301  
##  Residual             0.3472   0.5893  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)                           -1.16381    0.24347  -4.780 1.75e-06 ***
## tree_nodes                             0.23014    0.02589   8.889  < 2e-16 ***
## tree_health_simplifiedozone            0.20237    0.09565   2.116  0.03436 *  
## tree_health_simplifiedozone_and_other  0.21658    0.07918   2.735  0.00623 ** 
## tree_expositionexposed                 0.11549    0.06838   1.689  0.09123 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____
## tree_nodes  -0.916                     
## tr_hlth_smp -0.076 -0.059              
## tr_hlth_s__ -0.189 -0.009  0.476       
## tr_xpstnxps -0.034 -0.063 -0.092  0.072

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_4_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested and tree exposition data
glm_5_coh<-glmer(tree_heigth ~ tree_nodes + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_5_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes + (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1537.7   1553.7   -764.8   1529.7      398 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5213 -0.6705 -0.1426  0.4839  4.9989 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.1129   0.3360  
##  Residual             0.3535   0.5946  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept) -1.03026    0.24273  -4.245 2.19e-05 ***
## tree_nodes   0.23528    0.02619   8.984  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## tree_nodes -0.939

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_5_coh)

########################
glm_6_coh<-glmer(tree_heigth ~ tree_nodes*tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.149233 (tol = 0.002, component 1)

summary(glm_6_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes * tree_health_simplified + tree_exposition +  
##     reforested + (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1529.3   1569.2   -754.6   1509.3      392 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5718 -0.6976 -0.0957  0.5028  4.9427 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.09461  0.3076  
##  Residual             0.33090  0.5752  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                                   Estimate Std. Error t value
## (Intercept)                                      -1.134342   0.384100  -2.953
## tree_nodes                                        0.237738   0.042440   5.602
## tree_health_simplifiedozone                       0.272458   0.619126   0.440
## tree_health_simplifiedozone_and_other             0.399009   0.464316   0.859
## tree_expositionexposed                            0.113422   0.067966   1.669
## reforestedyes                                    -0.323423   0.104602  -3.092
## tree_nodes:tree_health_simplifiedozone           -0.005865   0.070827  -0.083
## tree_nodes:tree_health_simplifiedozone_and_other -0.021693   0.052639  -0.412
##                                                  Pr(>|z|)    
## (Intercept)                                       0.00314 ** 
## tree_nodes                                       2.12e-08 ***
## tree_health_simplifiedozone                       0.65989    
## tree_health_simplifiedozone_and_other             0.39015    
## tree_expositionexposed                            0.09515 .  
## reforestedyes                                     0.00199 ** 
## tree_nodes:tree_health_simplifiedozone            0.93400    
## tree_nodes:tree_health_simplifiedozone_and_other  0.68026    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____ tr_xps rfrstd tr_:__
## tree_nodes  -0.965                                          
## tr_hlth_smp -0.514  0.516                                   
## tr_hlth_s__ -0.749  0.753  0.433                            
## tr_xpstnxps -0.070  0.011  0.032  0.069                     
## reforestdys -0.175  0.099  0.125  0.045  0.021              
## tr_nds:tr__  0.519 -0.532 -0.988 -0.432 -0.047 -0.138       
## tr_nds:____  0.741 -0.767 -0.432 -0.986 -0.058 -0.043  0.443
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.149233 (tol = 0.002, component 1)

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_6_coh)

########################
glm_7_coh<-glmer(tree_heigth ~ tree_nodes*tree_health_simplified + tree_exposition + reforested + (1|plot), data = condition_model_all_C, family = Gamma(link = "log"))

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.149233 (tol = 0.002, component 1)

summary(glm_7_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes * tree_health_simplified + tree_exposition +  
##     reforested + (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1529.3   1569.2   -754.6   1509.3      392 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5718 -0.6976 -0.0957  0.5028  4.9427 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.09461  0.3076  
##  Residual             0.33090  0.5752  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                                   Estimate Std. Error t value
## (Intercept)                                      -1.134342   0.384100  -2.953
## tree_nodes                                        0.237738   0.042440   5.602
## tree_health_simplifiedozone                       0.272458   0.619126   0.440
## tree_health_simplifiedozone_and_other             0.399009   0.464316   0.859
## tree_expositionexposed                            0.113422   0.067966   1.669
## reforestedyes                                    -0.323423   0.104602  -3.092
## tree_nodes:tree_health_simplifiedozone           -0.005865   0.070827  -0.083
## tree_nodes:tree_health_simplifiedozone_and_other -0.021693   0.052639  -0.412
##                                                  Pr(>|z|)    
## (Intercept)                                       0.00314 ** 
## tree_nodes                                       2.12e-08 ***
## tree_health_simplifiedozone                       0.65989    
## tree_health_simplifiedozone_and_other             0.39015    
## tree_expositionexposed                            0.09515 .  
## reforestedyes                                     0.00199 ** 
## tree_nodes:tree_health_simplifiedozone            0.93400    
## tree_nodes:tree_health_simplifiedozone_and_other  0.68026    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____ tr_xps rfrstd tr_:__
## tree_nodes  -0.965                                          
## tr_hlth_smp -0.514  0.516                                   
## tr_hlth_s__ -0.749  0.753  0.433                            
## tr_xpstnxps -0.070  0.011  0.032  0.069                     
## reforestdys -0.175  0.099  0.125  0.045  0.021              
## tr_nds:tr__  0.519 -0.532 -0.988 -0.432 -0.047 -0.138       
## tr_nds:____  0.741 -0.767 -0.432 -0.986 -0.058 -0.043  0.443
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.149233 (tol = 0.002, component 1)

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_7_coh)

########################
# Model all independent variables WITH health status interaction WITHOUT reforested data
glm_8_coh<-glm(tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition + reforested, data = condition_model_all_C, family = Gamma(link = "log"))
summary(glm_8_coh)

## 
## Call:
## glm(formula = tree_heigth ~ tree_nodes + tree_health_simplified + 
##     tree_exposition + reforested, family = Gamma(link = "log"), 
##     data = condition_model_all_C)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7444  -0.6398  -0.1405   0.3019   1.7701  
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                           -0.63394    0.24337  -2.605  0.00954 ** 
## tree_nodes                             0.19493    0.02629   7.416 7.41e-13 ***
## tree_health_simplifiedozone            0.23272    0.10299   2.260  0.02438 *  
## tree_health_simplifiedozone_and_other  0.18895    0.07914   2.388  0.01742 *  
## tree_expositionexposed                 0.13374    0.07557   1.770  0.07751 .  
## reforestedyes                         -0.32034    0.07950  -4.029 6.71e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.4626594)
## 
##     Null deviance: 235.07  on 401  degrees of freedom
## Residual deviance: 185.21  on 396  degrees of freedom
## AIC: 1615.5
## 
## Number of Fisher Scoring iterations: 6

#### PLOT MODEL
par(mfrow =c(2,2))
plot(glm_8_coh)

# Choose better model
AIC(glm_1_coh,glm_2_coh,glm_3_coh, glm_4_coh, glm_5_coh, glm_6_coh, glm_7_coh, glm_8_coh)

######
library(sjmisc)
library(lmerTest)
library(bruceR)
summary(glm_2_coh)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: tree_heigth ~ tree_nodes + tree_health_simplified + tree_exposition +  
##     reforested + (1 | plot)
##    Data: condition_model_all_C
## 
##      AIC      BIC   logLik deviance df.resid 
##   1525.5   1557.5   -754.7   1509.5      394 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5686 -0.7133 -0.0860  0.4988  4.9872 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot     (Intercept) 0.09391  0.3064  
##  Residual             0.33133  0.5756  
## Number of obs: 402, groups:  plot, 45
## 
## Fixed effects:
##                                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)                           -1.02639    0.24521  -4.186 2.84e-05 ***
## tree_nodes                             0.22537    0.02568   8.776  < 2e-16 ***
## tree_health_simplifiedozone            0.22387    0.09537   2.348   0.0189 *  
## tree_health_simplifiedozone_and_other  0.21017    0.07834   2.683   0.0073 ** 
## tree_expositionexposed                 0.11226    0.06787   1.654   0.0981 .  
## reforestedyes                         -0.32335    0.10358  -3.122   0.0018 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tr_nds tr_hl_ tr____ tr_xps
## tree_nodes  -0.913                            
## tr_hlth_smp -0.054 -0.071                     
## tr_hlth_s__ -0.195 -0.006  0.468              
## tr_xpstnxps -0.035 -0.063 -0.093  0.071       
## reforestdys -0.184  0.064 -0.073  0.021  0.015

r2(glm_2_coh)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.429
##      Marginal R2: 0.267

model_performance(glm_2_coh)

## Report glmer

Figure S11

print_table(glm_2_coh)

## ────────────────────────────────────────────────────────────────────────
##                                        Estimate    S.E.      t     p    
## ────────────────────────────────────────────────────────────────────────
## (Intercept)                              -1.026 (0.245) -4.186 <.001 ***
## tree_nodes                                0.225 (0.026)  8.776 <.001 ***
## tree_health_simplifiedozone               0.224 (0.095)  2.348  .019 *  
## tree_health_simplifiedozone_and_other     0.210 (0.078)  2.683  .007 ** 
## tree_expositionexposed                    0.112 (0.068)  1.654  .098 .  
## reforestedyes                            -0.323 (0.104) -3.122  .002 ** 
## ────────────────────────────────────────────────────────────────────────

model_summary(glm_2_coh)

## 
## Model Summary
## 
## ──────────────────────────────────────────────────────
##                                        (1) tree_heigth
## ──────────────────────────────────────────────────────
## (Intercept)                              -1.026 ***   
##                                          (0.245)      
## tree_nodes                                0.225 ***   
##                                          (0.026)      
## tree_health_simplifiedozone               0.224 *     
##                                          (0.095)      
## tree_health_simplifiedozone_and_other     0.210 **    
##                                          (0.078)      
## tree_expositionexposed                    0.112       
##                                          (0.068)      
## reforestedyes                            -0.323 **    
##                                          (0.104)      
## ──────────────────────────────────────────────────────
## Marginal R^2                              0.267       
## Conditional R^2                           0.429       
## AIC                                    1525.487       
## BIC                                    1557.459       
## Num. obs.                               402           
## Num. groups: plot                        45           
## Var: plot (Intercept)                     0.094       
## Var: Residual                             0.331       
## ──────────────────────────────────────────────────────
## Note. * p < .05, ** p < .01, *** p < .001.
## 
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                    Term  VIF    VIF 95% CI Increased SE Tolerance
##              tree_nodes 1.02 [1.00, 10.75]         1.01      0.98
##  tree_health_simplified 1.04 [1.00,  1.49]         1.02      0.96
##         tree_exposition 1.03 [1.00,  1.77]         1.02      0.97
##              reforested 1.01 [1.00, 29.99]         1.01      0.99
##  Tolerance 95% CI
##      [0.09, 1.00]
##      [0.67, 1.00]
##      [0.56, 1.00]
##      [0.03, 1.00]

Fig 5a: Plot model aLL data with cohort

#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all_C,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5a<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "inverse")), aes(color= tree_health_simplified, fullrange =T))+
  scale_color_manual(values=  my_cols, labels= spanish_labels,
                    name= "Estado de salud")+ 
  scale_shape_manual(name= "Estado de salud", spanish_labels, values = c(15,16,17)) + 
  labs(x="Edad del árbol (años)", y= "Altura de los árboles")+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("a)")

## Warning in geom_smooth(method = "glm", method.args = list(family = Gamma(link =
## "inverse")), : Ignoring unknown aesthetics: fullrange

fig5a

Fig 5b: Plot model ORIGIN data with cohort

#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

fig5b<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
  ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified, linetype = reforested), fullrange =T)+
  scale_color_manual(values=  my_cols, labels= spanish_labels,
                    name= "Estado de salud")+ 
  scale_shape_manual(name="Estado de salud", labels=spanish_labels, values = c(15,16,17)) + 
  scale_linetype_manual(values= c("solid", "dashed"),labels = c("Regeneración\n natural", "Reforestados"), name= "Origen")+
  labs(x="Edad del árbol (años)", y= "Altura de los árboles")+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("b)")
fig5b

Fig 5c: model EXPOSITION data with cohort

#
# https://rpubs.com/dgolicher/ggplots2
# https://stackoverflow.com/questions/46302410/how-to-specify-link-function-in-ggplot-glm-graph

g0<-ggplot(condition_model_all_C,aes(x=tree_nodes,y=tree_heigth))

# method = "glm", method.args = list(family = Gamma(link = "log"))
fig5c<-g0+geom_point(aes(colour=tree_health_simplified,
                         shape=tree_health_simplified), alpha=0.5, size= 2.5)+
    ylim(0,15)+
  geom_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), aes(color= tree_health_simplified,linetype = tree_exposition), fullrange =T)+
  scale_color_manual(values=  my_cols, labels= spanish_labels,
                    name= "Altura de los árboles")+
  scale_shape_manual(name= "Altura de los árboles", labels=spanish_labels, values = c(15,16,17)) + 
  scale_linetype_manual(values= c("solid", "dashed"),labels = c("Con nodriza", "Expuestos"), name= "Exposición de los árboles",) +
  labs(x="Edad del árbol (años)", y= "Altura de los árboles")+
  #stat_smooth(method = "glm", method.args = list(family = Gamma(link = "log")), se = TRUE)+
  theme_bw()+ 
  theme(legend.key.size = unit(0.7, 'cm'),
        legend.key.width= unit(0.7, 'cm'),
        legend.title = element_text(size=10),
        legend.text = element_text(size=10))+
  theme(text = element_text(size = 20))+
  theme(legend.title.align = 0.5)+
  ggtitle("c)")
fig5c

Multiplot Fig 5

multiplot(fig5a, fig5b, fig5c, cols=1)

Figure S8

p_d <- ggplot(parcelas_long, aes(x=plot, y=n_trees, fill=tree_health_simplified)) +
  geom_bar_pattern(stat="identity",
                   aes(pattern = tree_health_simplified),
                   color = "black",
                   pattern_fill = "black",
                   pattern_density = 0.1,
                   pattern_spacing = 0.015,
                   pattern_key_scale_factor = 0.6) +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") + 
  scale_pattern_manual(name= "Estado de salud", breaks = desired_order,
                       values=c ("none",  "none", "none", "none", 
                                 "crosshatch", "none", "stripe", "circle", "none", "none"),
                       labels = spanish_labels)
  

p_d <- p_d + theme_bw() +
  labs(x="Parcelas", y= "Número de árboles") +
  theme(text = element_text(size = 20)) +
  ggtitle("")+
  coord_flip()
p_d 

Figure S12

# plot pies in map
p_satmap <-  ggmap(sat_map) 
p_satmap +geom_scatterpie(data=parcelas_tidy,
                aes(x=X_coordinates_longitude,
                    y=X_coordinates_latitude,
                    group=plot),
                pie_scale = 1.5,
                cols=desired_order,
                color=NA,
                alpha=1)  +
  scale_fill_manual(values= my_cols, breaks = desired_order,
                    labels= spanish_labels,
                    name= "Estado de salud") +
  theme(text = element_text(size = 20))

Figure S13

# Create new variable with percentage of ozone damage
parcelas_tidy<-parcelas_tidy %>% rowwise() %>% 
                     mutate(., 
                      total=sum(healthy,ozone,ozone_and_other,
                          drougth, acid_rain, other,
                          others_combined, dead, fungi,
                          # insect, 
                          worm)) %>%
                    mutate(perc.ozone= sum(ozone, ozone_and_other)/total)

#Statistical
leveneTest( parcelas_tidy$X_coordinates_altitude, parcelas_tidy$perc.ozone)

## Warning in leveneTest.default(parcelas_tidy$X_coordinates_altitude,
## parcelas_tidy$perc.ozone): parcelas_tidy$perc.ozone coerced to factor.

## Warning in anova.lm(lm(resp ~ group)): ANOVA F-tests on an essentially perfect
## fit are unreliable

#plot
p <- ggplot(parcelas_tidy) +
     geom_point(aes(x=X_coordinates_altitude,
             y=perc.ozone))


p<- ggplot(parcelas_tidy, aes(X_coordinates_altitude, perc.ozone ))+
  geom_point(color= "grey50", size = 3, alpha = 0.6)

p<-p + 
  stat_smooth(color = "skyblue", formula = y ~ x,fill = "skyblue", method = "lm") +
  stat_poly_eq(
    aes(label = paste(..eq.label.., ..adj.rr.label.., sep = '~~~~')),
    formula = y ~ x,  parse = TRUE,
      size = 5, # Tamaño de fuente de la fórmula
             label.x = 0.1, #location, la proporción entre 0-1
      label.y = 0.95)+
  labs(x="Altitud de la parcela", y= "Percentaje de daño de los árboles dentro de las parcelas")+
  theme_bw() +
   theme(plot.title = element_text(lineheight=1.1, face="bold")) +
  theme(text = element_text(size = 20), 
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        axis.text=element_text(size=9))

p

## Warning: The dot-dot notation (`..eq.label..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(eq.label)` instead.

ggsave("../outputs/FigS3.png",
       dpi = 600)

RESULTS Mean and sd of trees collected by each ranger during Participatory monitoring

# No se incluye a M (37) y PP (6), ya que eran colectados por diferentes miembros del team
trees<- c(135,133,130,174,150,139,165,166,133,137,125,105)
rangers<- c("A","B","C","D","E","F","G","H", "I", "J", "K", "L")
mean(trees)

## [1] 141

sd(trees)

## [1] 19.60519

Percentage damage in plots

#Percentage damage plots

PDP<-data.frame(parcelas_tidy$plot, parcelas_tidy$perc.ozone)

PDP<-PDP[order(PDP$parcelas_tidy.plot),]
range(PDP$parcelas_tidy.perc.ozone)

## [1] 0.00000 0.84375

head(PDP[order(PDP$parcelas_tidy.perc.ozone),])

Percentage per each health status

# Create new variable with porcentage of ozone damage

status_HS<- c("healthy","ozone","ozone_and_other",
                          "drougth", "acid_rain", "other",
                          "others_combined", "dead", "fungi",
                          "worm")
count_HS<-c(sum(parcelas_tidy$healthy), sum(parcelas_tidy$ozone),
            sum(parcelas_tidy$ozone_and_other),sum(parcelas_tidy$drougth),
            sum(parcelas_tidy$acid_rain),sum(parcelas_tidy$other),
            sum(parcelas_tidy$others_combined), sum(parcelas_tidy$dead),
            sum(parcelas_tidy$fungi), sum(parcelas_tidy$worm))
sum(count_HS)

## [1] 1765

# Tree number to each health status
percentageH<-data.frame(status_HS, count_HS)
percentageH

# Percentage of ozone and ozone + other trees
(sum(c(parcelas_tidy$ozone, parcelas_tidy$ozone_and_other))*100)/sum(count_HS)

## [1] 35.35411

# Percentage of healthy trees
(sum(parcelas_tidy$healthy)*100)/sum(count_HS)

## [1] 27.70538

# Percentage of others damages in trees with dead trees
(sum(c(parcelas_tidy$drougth, parcelas_tidy$acid_rain , parcelas_tidy$other, parcelas_tidy$others_combined,parcelas_tidy$dead, parcelas_tidy$fungi, parcelas_tidy$worm))*100)/sum(count_HS)

## [1] 36.94051

# Filtrar por categoría de daño
condition_HOz<-muestreo_tidy%>%
  filter(tree_health_simplified == "healthy" | tree_health_simplified == "ozone")

condition_HOz$tree_health_simplified <- as.factor(condition_HOz$tree_health_simplified)


# Damages by other NO ozone NO DEAD
levels(as.factor(muestreo_tidy$tree_health_simplified))

##  [1] "acid_rain"       "dead"            "drougth"         "fungi"          
##  [5] "healthy"         "other"           "others_combined" "ozone"          
##  [9] "ozone_and_other" "worm"

no_ozone<- dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "acid_rain"| tree_health_simplified == "fungi" | tree_health_simplified == "others_combined" | tree_health_simplified == "worm" | tree_health_simplified == "drougth"| tree_health_simplified == "other")

(609*100)/1765

## [1] 34.50425

# Percentage of DEAD trees
levels(as.factor(muestreo_tidy$tree_health_simplified))

##  [1] "acid_rain"       "dead"            "drougth"         "fungi"          
##  [5] "healthy"         "other"           "others_combined" "ozone"          
##  [9] "ozone_and_other" "worm"

dead<- dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "dead")
count(dead)

(43*100)/1765

## [1] 2.436261

sum(489,125, 499,1,9,70, 3,271,255, 43)

## [1] 1765

Chi test

count(dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "ozone"))

count(dplyr::select(muestreo_tidy, contains(c("tree_health_simplified", "reforested", "tree_exposition"))) %>%
  filter(tree_health_simplified == "ozone_and_other"))

OBS_ozo<-c(125,499)
ESP_ozo<-c(0.5, 0.5)
chisq.test(x = OBS_ozo,
           p = ESP_ozo)

## 
##  Chi-squared test for given probabilities
## 
## data:  OBS_ozo
## X-squared = 224.16, df = 1, p-value < 2.2e-16

#224.1603
#p-value < 2.2e-16
# p<0.001